home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / os2 / adaptor.zip / ADAPT.ZIP / adaptor / examples / f77 / lumm.f < prev    next >
Text File  |  1993-06-28  |  3KB  |  141 lines
  1.       PROGRAM LUMM
  2. c     IMPLICIT NONE
  3.       INTEGER N
  4.       REAL A(:,:), lu(:,:), b(:,:)
  5.       REAL COL (:)
  6. cmf$  layout col (:serial)
  7.       INTEGER I,J, K
  8.       real maxdiff, eps
  9.       real rand
  10.       real x
  11.  
  12.       print *, 'What is the size of the Matrix : '
  13.       read *, N
  14.  
  15.       allocate (a(1:N,1:N), b(1:N,1:N), lu(1:n,1:n), col(n))
  16.  
  17.       print *, 'Arrays are allocated'
  18.  
  19.       call cmf_random (A)
  20.  
  21. !HPF$ INDEPENDENT, LOCAL_ACCESS
  22.       DO J=1,N 
  23. !HPF$ INDEPENDENT, LOCAL_ACCESS
  24.           DO I=1,N
  25.              lu (i,j) = a(i,j)
  26.           END DO
  27.       END DO 
  28.  
  29.       PRINT *,'PROGRAM STARTS'
  30.       DO K=1,N-1
  31.          x = lu(k,k)
  32.          do j = k+1, n
  33.             lu(j,K) = lu(j,K)/x
  34.          end do
  35.          col = lu(1:n,k)
  36. !HPF$ INDEPENDENT, LOCAL_ACCESS
  37.          do  j = k+1, n
  38. !HPF$ INDEPENDENT, LOCAL_ACCESS
  39.             do i = k+1, n
  40.               lu (i,j) = lu(i,j) - col(i) * lu(k,j)
  41.             end do
  42.          end do
  43.       END DO
  44.  
  45.       PRINT *, 'LU ready'
  46.       PRINT *, 'now finds Y with L * Y = E '
  47.  
  48. !HPF$ INDEPENDENT, LOCAL_ACCESS
  49.       DO J=1,N
  50. !HPF$ INDEPENDENT, LOCAL_ACCESS
  51.           DO I=1,N
  52.               if (I .eq. J) then
  53.                  B(I,J) = 1
  54.                else
  55.                  B(I,J) = 0
  56.               end if 
  57.           END DO
  58.       END DO 
  59.  
  60.       do j=1, n
  61.          col = lu(1:n,j)
  62. !HPF$ INDEPENDENT, LOCAL_ACCESS
  63.          do k = 1,n
  64. !HPF$ INDEPENDENT, LOCAL_ACCESS
  65.             do i = j+1,n
  66.               b(i,k) = b(i,k) - col(i)*b(j,k)
  67.             end do
  68.          end do
  69.       end do
  70.  
  71.       PRINT *, 'now finds X with R * X = E '
  72.  
  73.       do j=n, 1, -1
  74.          col = lu(1:n,j)
  75. !HPF$ INDEPENDENT, LOCAL_ACCESS
  76.          do k = 1,n
  77.             b(j,k) = b(j,k) / col(j)
  78. !HPF$ INDEPENDENT, LOCAL_ACCESS
  79.             do i = 1, j-1
  80.               b(i,k) = b(i,k) - col(i) *b(j,k)
  81.             end do
  82.          end do
  83.       end do
  84.  
  85.       PRINT *, 'b is the inverse '
  86.  
  87. c     print *, 'a = ', a
  88. c     print *, 'b = ', b
  89.  
  90.       PRINT *, 'lu = a * b '
  91.  
  92. c     lu = 0
  93. !HPF$ INDEPENDENT, LOCAL_ACCESS
  94.       DO J=1,N
  95. !HPF$ INDEPENDENT, LOCAL_ACCESS
  96.           DO I=1,N
  97.              lu (i,j) = 0
  98.           END DO
  99.       END DO 
  100.  
  101.       do k=1,n
  102.          col = a(1:n,k)
  103. !HPF$ INDEPENDENT, LOCAL_ACCESS
  104.          do j=1,n
  105. !HPF$ INDEPENDENT, LOCAL_ACCESS
  106.            do i=1,n
  107.              lu(i,j) = lu(i,j) + col(i) * b(k,j)
  108.            end do
  109.          end do
  110.       end do
  111.           
  112.       PRINT *, 'matmul correct, now check = E'
  113.  
  114. c     print *, 'E = ', lu
  115.  
  116. !HPF$ INDEPENDENT, LOCAL_ACCESS
  117.       do i=1,N
  118. !HPF$ INDEPENDENT, LOCAL_ACCESS
  119.          do j=1,N
  120.             if (i .eq. j) then
  121.                a(i,j) = lu(i,j) - 1
  122.              else
  123.                a(i,j) = lu(i,j)
  124.             end if
  125.             a(i,j) = abs (a(i,j))
  126.          end do
  127.       end do
  128.  
  129.       maxdiff = 0.0
  130. !HPF$ INDEPENDENT, LOCAL_ACCESS
  131.       do i=1,N
  132. !HPF$ INDEPENDENT, LOCAL_ACCESS
  133.          do j=1,N
  134.            REDUCE (maxval, maxdiff, a(i,j))
  135.          end do
  136.       end do
  137.  
  138.       deallocate (col, lu, b, a)
  139.       print * , 'maximal eps = ', maxdiff
  140.       END
  141.